program VMD !  works with df   based on code drivenDf Runge K.f95 which is in frame of source (phase velocity=0)   

  use fft_wrap, only   : fft_init, fft_forward, fft_backward

  use dimensions, only : jx, jz, jv, jDirection, &   ! domain size in r-space
                         nkx, nkz, np,           &   ! limits      in k-space
                         dx,  dz,  dv,           &   ! resolution  in r-space
                         dkx, dkz, dp                ! resolution  in k-space
         

  use parameters, only : phi0, kz, vPhase,  Dxxxx, Dvvvv, &            ! real,     phys
                         ny,  iWaveFrame, iFixPump, iModx, iModz, &    ! integer,  phys
                         ModAmp,  Lx, SpeckleWidth, &                  ! real,     phys
                         restart, nt, &                                ! integer,  num
                         Vmax, dt, KxFilter, &                         ! real,     num 
                         iself, iskip, idfMod, &                       ! integer,  num 
                         iE0Modulation, iBcPerpPeriodic, iFourier, &   ! integer,  num
                         pi, iunit, isign                              ! constants


  implicit none



  !-- variables  ----------------------- 

  integer    :: ix, iz, iv, iDirection  ! counters in r-space  
  integer    :: ikx, ikz, ip            ! counters in k-space
  integer    :: it, itMin, itMax        ! iteration count
 
  real       :: OmegaBounce

  real       :: t, u, v0

  real       :: Courant_in_x,  Courant_in_z   

  complex    :: ExpKVdt, phi0T

   

  !-- arrays -------------------------


  real       :: dfu(jz, jv, jx, jDirection)
  real       :: uxMat(jx, jz, jDirection)
  real       :: uzMat(jx, jz, jDirection)

  real       :: Ez2D(jx,jz), Ex2D(jx,jz)
  real       :: E0z(jx,jz), Ezinternal(jx,jz)
  real       :: Density2Dreal(jx,jz)

  complex    :: phi02DFT(jx,jz)
  complex    :: phi2DFT(jx,jz)
  complex    :: CEz2D(jx,jz)

  complex    :: EzKx(jx,jz)
  complex    :: ExpKv(jz,jv)
  complex    :: ExpKx(jx,jDirection)


  !-- vectors -------------------------

  real       :: xvec(jx), zvec(jz), vvec(jv)      ! coordinate vectors in r-space
  real       :: kXvec(jx), kZvec(jz), pvec(jv)    ! coordinate vectors in k-space

  real       :: ExFilter(jx)

  real       :: avEx2(jx)

  real       :: f0vec(jv)
   
  complex    :: xMod(jx), uPerpMod(jx)
  complex    :: density1DFT(jz)
  complex    :: E0C(jz), Phi0C(jz)

  complex    :: fvc0T(jv)
  complex    :: ExpKz(jz)


  data       it, itMin / 0, 1 /
  data       t, u /0.0, 1.0/

 
!-------------------------------------------------------------------------------------


      call fft_init()

      call vmd_init(u, v0, xvec, zvec, vvec, kXvec, kZvec, pvec)       ! outputs


      !-- reading restart file or create new initial conditions --
      
      call vmd_ic(xvec, zvec, vvec, kZvec, &                           ! inputs
                  f0vec, dfu, uxMat, E0C, phi0C, phi0T, fvc0T, t)      ! outputs


      !-- creating filter and bunch of "exp" vectors ---

      call vmd_init_ex(kXvec, kZvec, pvec, vvec, u, v0, &              ! inputs
                       ExFilter, ExpKVdt, ExpKz, ExpKv, ExpKx)         ! outputs

      
      
      !-- creating two modulation vectors --

      call vmd_init_modulation(xvec, kXvec,    &                       ! inputs 
                               xMod, uPerpMod)                         ! outputs


      !-- creating initial density, its FFT and modulated E field --

      call vmd_init_density(zvec, ExFilter, xMod, uPerpMod, E0C,  &    ! inputs
                            E0z,  phi02DFT)                            ! outputs


      !-- compute intitial electric field --
      
      call GetEfield(dfu, kXvec, kZvec, ExFilter, phi02DFT, &          ! inputs 
                     Ez2D, Ex2D, OmegaBounce, Density2Dreal, &         ! outputs
                     density1DFT, phi2DFT, &                           ! outputs
                     CEz2D, EzKx)                                      ! outputs



      !-- diagnostics --

      call write_csi()

      !-- write initial electric field --

      call write2D_vs_xz("E0z_initial_vs_xz.dat", E0z, xvec, zvec)
      call write2D_vs_xz("Ex_initial_vs_xz.dat", Ex2D, xvec, zvec)
      call write2D_vs_xz("Ez_initial_vs_xz.dat", Ez2D, xvec, zvec)


      !-- write headers for files recording time evolution  --


      write(*,*)
      if(iself==1)then
         write(*,*)'it/nt            t         Phi Internal   OmegaBounce     potEx         potEz'
      else
         write(*,*)'it/nt          kt        real(Epsilon)   Total Energy     Kinetic Energy'
      endif


      open(122,file='Ez_vs_kxtime.dat')
      write(122,*) '% t\kx'
      write(122,*) 'NA', ( kxVec(ikx), ikx = 1,nkx/2  - nkx/6 -1 )


      open(3,file="vstime.txt")
      if(iself==1)then
         write(3,*) &
      '%t   PhiInternal  PhiToTal -Phase(phi)  KineticEnergy Ex^2energy  Ez^2energy  OmegaBounce Courant_in_z'
      else
         write(3,*) &
      '%t   1/|Epsilon|  Re(Epsilon)  Im(Epsilon)  Re(Csi) Im(Csi) KineticEnergy  Phi0Energy   TotalEnergy   f^2 ntrapped'
      endif
      close(3)


!-------------------------------------------------------------------------------------
! BEGIN EVOLUTION
!-------------------------------------------------------------------------------------


 
do it=itMin,itMIn + nt -1  ! BEGIN EVOLUTION



   call evolve(   dfu, uxMat, uzMat, Ex2D, Ez2D,  E0C, f0vec, &    ! input/output
                  vvec, pvec, zvec, kXvec, kZvec, &                ! input
                  ExFilter, phi02DFT,  fvc0T, &                    ! input
                  xMod, uPerpMod, &                                ! input
                  ExpKv, ExpKz, ExpKVdt, &                         ! input
                  E0z, Density2Dreal, &                            ! output
                  phi2DFT,  &                                      ! output
                  CEz2D, EzKx,  density1DFT, &                     ! output
                  OmegaBounce )                                    ! output

   phi0T =    ExpKVdt * ExpKVdt * phi0T
   phi0C(:) = ExpKVdt * ExpKVdt * phi0C(:)


   t=t+dt

      
   Courant_in_z = maxval(abs(uzMat))*dt/dz
   Courant_in_x = maxval(abs(uxMat))*dt/dx

   !-- diagnostics --

   if(  ((it/iskip)*iskip.eq.it) .or. (it.eq.itMIn + nt -1) .or. (it.eq.itMin)  ) then      

      Ezinternal = Ez2D - E0z


      call diagnostics( t, float(it-itMin)/nt, Courant_in_z, &   ! inputs
                        OmegaBounce, SpeckleWidth,  &            ! inputs
                        dfu, vvec, f0vec, &                      ! inputs
                        Ezinternal, Ez2D, Ex2D, density1DFT)     ! inputs

      do ikx=1,nkx/2  - nkx/6 -1
         avEx2(ikx) = 0
         do iz = 1, jz
            avEx2(ikx) = avEx2(ikx) + abs( EzKx(ikx,iz) )**2
         end do
      end do
      write(122,*) t, ( avEx2(ikx), ikx = 1, nkx/2  - nkx/6 -1 )

   end if


      
   if(  Courant_in_z > 1.0 ) then
      write(*,*)   'Courant for z advection =', Courant_in_z, '  is > 1, stop evolution'
      go to 314
   endif
            
   if(  Courant_in_x > 1.0 ) then
      write(*,*)   'Courant for x advection =', Courant_in_x, '  is > 1, stop evolution'
      go to 314
   endif

end do      ! END TEMPORAL LOOP
 
314 continue

     close(122)  ! 'Ez_vs_kxtime.dat'

     open(5, file='output', position='append')
     write(5,*) 'Dvvvv should be > E_max / p_max^3 = ', maxval(ez2d) / pvec(np/2)**3
     close(5) 

     itMax = it
 
     write(*,*) 'end temporal evolution'  


!-------------------------------------------------------------------------------------
! END TEMPORAL LOOP
!-------------------------------------------------------------------------------------


     Ezinternal = Ez2D - E0z


     !-- writing 2D data as function of (x,z) --


     call write2D_vs_xz("E0z_final_vs_xz.dat",   E0z, xvec, zvec)

     call write2D_vs_xz("Ez_internal_vs_xz.dat", Ezinternal, xvec, zvec)
     call write2D_vs_xz("Ez_total_vs_xz.dat",    Ez2D,       xvec, zvec)
     call write2D_vs_xz("Ex_vs_xz.dat",          Ex2D,       xvec, zvec)

     call write2D_vs_xz("uxPlus_vs_xz.dat",  uxMat(:,:,1)-u,  xvec, zvec)
     call write2D_vs_xz("uxMinus_vs_xz.dat", uxMat(:,:,2)+u,  xvec, zvec)
     call write2D_vs_xz("uzPlus_vs_xz.dat",  uzMat(:,:,1)-v0, xvec, zvec)
     call write2D_vs_xz("uzMinus_vs_xz.dat", uzMat(:,:,2)-v0, xvec, zvec)

     call write2D_vs_xz("dnTotal_vs_xz.dat", Density2Dreal,   xvec, zvec)


     do iDirection=1,2    ! caution: this overwrites "density2Dreal" 
     
       do ix=1,jx
         do iz=1,jz
            density2Dreal(ix,iz) = sum(dfu(iz,:,ix,iDirection)*dv) 
         end do
       enddo
     
       if(iDirection==1)then
          call write2D_vs_xz("dnPlus_vs_xz.dat",  Density2Dreal,   xvec, zvec)
       else
          call write2D_vs_xz("dnMinus_vs_xz.dat", Density2Dreal,   xvec, zvec)
       endif
    
     enddo 

 
    !-- writing 2D distribution function, f(z,v) ---  

    call write_fdist(dfu, zvec, vvec, OmegaBounce)         ! inputs



    !-- writing some integral quantities ---  


    call write_means(dfu, kZvec, pvec, fvc0T)              ! inputs


    !-- write some specrtral spectrum quantities --


    call write_spectral(uzMat, uxMat, kXvec, kZvec, u,             &  ! inputs 
                        density1DFT, EzKx, CEz2D, phi2DFT)            ! inputs
 

    !-- write restart file --
  
    write(*,*)'phi0=',phi0,'    kz=',kz,' vPhase=',vPhase,'   Lx=', Lx,'     SpeckleWidth=',SpeckleWidth
    write(*,*)'begin restart file write' 

    open(4,file="restart")
    write(4,*) jz,jv,t,ny,phi0,kz,vPhase,ny,Dxxxx,Dvvvv,iWaveFrame, iFixPump, iModx,iModz, ModAmp, Lx, SpeckleWidth
    write(4,*) E0C, phi0C, phi0T      
    write(4,*) uxMat 
    write(4,*) dfu
    close(4)

    write(*,*)'restart file written'; write(*,*)
      

501   format(  6400(e12.6,2x)  )
       
   stop;

end


!-------------------------------------------------------------------------------------
! subroutine deAlias
!-------------------------------------------------------------------------------------

      
subroutine deAlias(jz,et); implicit none;  integer i,jz;  complex et(jz)
if( jz<6 ) return
       do i=jz/2  - jz/6, jz/2 + 2 + jz/6;   et(i)=(0.0,0.0);  end do
return;  end



!-------------------------------------------------------------------------------------
